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Abstract. In this article, several aspects of the dynamics of a toy model for long- 
range Hamiltonian systems are tackled focusing on linearly unstable unmagnetized 
(i.e. force- free) cold equilibria states of the Hamiltonian Mean Field (HMF). For 
special cases, exact finite-A linear growth rates have been exhibited, including, in 
some spatially inhomogeneous case, finite- N corrections. A random matrix approach 
is then proposed to estimate the finite- A growth rate for some random initial states. 
Within the continuous, N — > oo, approach, the growth rates are finally derived without 
restricting to spatially homogeneous cases. All the numerical simulations show a very 
good agreement with the different theoretical predictions. Then, these linear results 
are used to discuss the large-time nonlinear evolution. A simple criterion is proposed 
to measure the ability of the system to undergo a violent relaxation that transports it 
in the vicinity of the equilibrium state within some linear e-folding times. 
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1. Introduction 

The Hamiltonian Mean Field (HMF) model has become a well-known toy model to 
address the intricate relationships between dynamics and statistical mechanics of long- 
range interacting systems. It is defined by the following Hamiltonian 

N 2 i N 

^ = ET + 2iv^ [1 - cos(ai -^ )] ' (1) 

1=1 i,i=l 

where N is the number of particles, and 6^ and Pi denote respectively the position 
and momentum of the i th particle. A useful collective quantity to introduce is the 
magnetization vector (M x ,M y ) with 

1 N 1 N 

M * = nJ2 cosdi and M v = n J2 sind * ( 2 ) 

1=1 1=1 
The average energy per particle U = H/N reads then 

N 2 i 

U = Y,h + \( l - M2 ), (3) 

i=l 

where M = M x 2 + M y 2 denotes the modulus of the magnetization vector. 

Recently, much interest has been devoted to the so-called quasi-stationary states 
(QSS) which are known to be responsible for the very slow convergence towards the 
statistical mechanics equilibrium predictions. Far from being difficult to generate, 
these QSS naturally emerge from waterbag initial distributions [HE]- It is also known 
that initial waterbag conditions in momenta, associated to zero or almost zero initial 
magnetization, induce the longest lasting QSS. However, when lowering towards zero 
the initial temperature of the particles, it is possible to exhibit waterbag momenta 
configurations with vanishing magnetization which converge exponentially towards the 
Boltzmann-Gibbs equilibrium. This calls for a linear theory approach. 

Linear stability of the HMF model about unmagnetized equilibrium states has 
been up to now mostly studied within the Vlasov framework [31 HO], which assumes 
in particular an infinite number of particles. Moreover, to the authors' knowledge, 
the linear stability of spatially inhomogeneous, unmagnetized, equilibria has not been 
considered yet, even in the continuous Vlasov approach. 

The motivation of the present study is then twofold: Firstly and mostly, one wishes 
to tackle the linear study of the unmagnetized cold HMF equilibria, within a finite- N, 
therefore exact, framework; secondly, the ensuing nonlinear dynamics is briefly addressed 
to show that the thermalization of cold unmagnetized HMF systems finely illustrates 
Lynden-Bell's concept of violent relaxation for long-range systems. 

In Section [2j we shall establish the finite- iV framework used for the linear stability 
derivation. In Section El we shall calculate the exact linear growth rates for two finite- 
iV equilibria, both of zero temperature and zero magnetization, and compare them 
to numerical simulations. Section H] is dedicated to a random matrix approach for 
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the calculation of non-deterministic symmetric initial states growth rates. Section 
tackles the general case, including the spatially inhomogeneous ones by introducing a 
fluid approach based on Vlasov framework. Section [6] ends this study by discussing 
the connections between the linear features just derived and the HMF thermalization 
properties. The dynamics of the cold unmagnetized HMF model is proposed as a 
paradigm of violent relaxation. 



2. Linear dynamics about cold unmagnetized finite- N equilibria 



The equations of motion can straightforwardly be written from Equation (Cp) as 

h = Pk 



Vke{l,...,N} 



N 



Pk 



Lj2sm(Oi-o k 



(4) 



N 



i=l 



Using Equation (J2J), the force acting on the particle k may be written as 

F k = M y cos k - M x sin k . (5) 

Let us consider unmagnetized finite-iV equilibria, namely stationary states of the 
equations of motion 01]), with M x = M y = 0. This amounts to have p\ = and 
angles 0% distributed in such a way that M x = M y = 0. Let us remark that, since the 
total momentum P = Yli=i Pi * s a constant of motion, the zero momentum equilibrium 
case considered here is just the cold (i.e. monokinetic) case for the special choice P = 0. 

Let us perform the linear stability of this system. In that purpose, we write 
9 k = 6* k + 59 k and p k = 5p k where the asterisk denotes the unperturbed solution. At 
first order in 59 %, the force SF k felt by the particle k verifies 

N N 



5F k 



1 

N 



cos 6* k ^2 cos O^SOi + sin e* k ^2 sin Ot^Oi 



i=l 



8=1 



This yields the following linear system 





S9 k 






On In 




' {80k} ' 




5pk 






A N 




{Spk} 



where In is the N x N identity matrix and A a, N x N matrix defined by 



M 4 = — cos( 



3)- 



(6) 



(7) 



(8) 



The stability depends on the eigenvalues {A^} of the Jacobian matrix in Equation ([7]). 
Let us name it J. The eigenvalue problem can now be reduced to the unique study of 
A using the following transformation 



J- XI 



2N 



-XIn On 
A I N 



In 
On 



1 T 

-XI n + \A 



(9) 



provided A 7^ 0. Hence, 

det( J-XIn) = det(-AIiv) det(-XI N - 



-A) = (-l) Ar det (A - X 2 I N ) -(10) 
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A,j = ~r COS ( — (i - j) ) . (12) 



In other words, writing xm the characteristic polynomial of M, one has 

Xj(X) = (-1) N Xa(X 2 ) (11) 

This factorization allows us to focus only on the matrix A, and deduce the eigenvalues 
of the higher-order matrix J by taking the square root of A ones. We shall now study 
the case of two particular finite- iV equilibria that advantageously simplify A. 

3. Exact finite-iV treatment for two special cold force-free equilibria 

3.1. The quiet start case 

Having in mind the computational plasma terminology, we define the so-called "quiet 
start" configuration as the equilibrium characterized by an equirepartition of the 
particles on the circle. It is here formally described by Wk, 91 = 27ik/N, up to 
some constant phase, and pi = 0. Using this definition in Equation ([8]), the matrix 
A immediately reduces to 

1 /2tt 
— cos — 
N \N 

We notice that we can rewrite the coefficients as 

A-i,3 = -^(i-i)modJV) (13) 

which makes evident that A is a circulant matrix. Using Equation ffTTj) . it can be easily 
shown that its eigenvalues {A^} are expressed as 

N 

\l = Y,A 3 e 2 ^l N . (14) 

Moreover, since A is real and symmetric, its eigenvalues are real, and one can identify 
the previous equation with its real part, yielding 

This is just 

A 2 fc = \ (4,i + 4,v-i) , (16) 

so that A has only one double non-zero eigenvalue equal to 1/2. Using Equation fflTT) 
finally yields the expected growth rate ■jqs coming from the spatially homogeneous 
Vlasov linear theory [3] for the cold waterbag as 

1QS = = 4, (17) 

with no finite-iV correction. In order to test the validity of this linear study, we 
performed numerical simulations based on a fourth-order symplectic integrator [3]. 
Starting from a quiet start configuration, every particle is moved by a uniformly 
randomized quantity e 2tt/N. As shown in Figured! the behaviour of the system 
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during the early times shows a very good agreement with the predicted exponential 
growth, and does not depend on the number of particles, which only changes the initial 
magnetization resulting from the perturbation. 




"0 2 4 6 8 10 

Time t 



Figure 1. Numerical plot of the magnetization with respect to time. As predicted by 
Equation (1171) . the slope does not depend on the number of particles N, and is in very 
good agreement with the theoretical value of l/\/2, plotted in thick dashed line. 



3.2. Bi-clustered quiet start 

The previous case involved a finite-iV analog of a homogeneous Vlasov force-free 
equilibrium. We can construct another finite- N equilibrium with zero magnetization 
by uniformly distributing N/2 particles in a cluster of size A8 centered on a given 
position, and by settling the positions of the N/2 remaining ones by rotating the first 
cluster by n with 

2 M e * N/2+k = e * h+ * (18) 

J 1 Pi = o 

Figure [2] shows an example of such a bi-cluster configuration. The same stability analysis 




Figure 2. Plot of some finite- N bi-clustered equilibrium where each particle has a 
zero momentum and faces its symmetric on the circle, providing a zero magnetization. 

as in the homogeneous quiet start case can be performed. However, one expects now 
the growth rate to depend on A9, the clusters' size. Equation f TTTj) is still valid, but 
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the matrix A is no longer circulant as it was in the simple quiet start case. A can be 
rewritten under the form 



A 



L -L 
-L L 



(19) 



where L is a N/2 x N/2 matrix with coefficients 



1 /4A0, A , s 

^ = iv cos (~ J' 

Therefore, the characteristic polynomial of A reads 

det (A - X 2 I N ) = (~2\ 2 ) N/2 det f L - y . (21) 



Or, equivalently, 



XA {\ 2 ) = {-2\ 2 ) N '\ L (T) (22) 



2 

This decomposition allows us to focus on the smaller matrix L. Unfortunately, L is not 
circulant either, but is a Toeplitz matrix. Indeed, one can write Ly = Luj\. A work 
performed by Treichler [6] showed that for a m x m Toeplitz matrix generated from the 
coefficients tk = cos(ku), the only two non-zero eigenvalues are 

/ \ 1 ( sin (moj) \ , , 

u±(m,u) = -[m± WM • 23 

Equation (123]) allows us to obtain the eigenvalues {u^} of the matrix L 

!~ ^ 4Ag > \ = I _ L sin(2A0) 
^ ~ N V± V 2 ' iV / 4 2iVsm (4A8/N) ' ^ ' 

Equations (jTTj) . (1221) and (I24I) give the growth rate as 



, I sin(2A#) 

7scqs - - y 2 + iV s in(4A^/iV)' (25) 

that is, in the large iV limit, 

y/1 + sinc(2Afl) 

Ibcqs = , (26) 

up to 0(N~ 2 ) terms. Figure ([3]) shows the comparison between the numerically 
computed growth rates and the theoretical prediction of Equation (|2"5|) . When A9 = 
7r/2, particles are uniformly distributed on the circle as in the previous equilibrium and 
the growth rate is exactly l/y2 with no finite- N correction. Otherwise, the growth 
rate depends on the number of particles but we checked that the difference with the 
asymptotic result (1261) is already very small for iV above 10 particles. 

We should note that in Reference [I], the authors diagonalize the matrix Ty = 
cos(#j — 9j) and give an analytic expression of its two non-zero eigenvalues. However, 
using their formula for a practical calculus requires the knowledge of an implicitely 
defined parameter ip, which analytical calculation, in a finite-iV context, is impossible, 
let alone some degenerate cases. 
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Figure 3. Plot of the growth rate as a function of A9 for slightly perturbed bicluster 
initial configurations. The circles correspond to exponential fits from the numerically 
integrated magnetization, and the solid curve corresponds to Equation (|25[) . Each 
numerical calculation has been performed with 1000 particles. 



4. Extension to symmetric random initial configurations 



In this section, we will only consider symmetric equilibria prepared to give M = 
in the following way: we distribute N/2 particles at random on a partition of a given 
length and put the N/2 remaining particles by shifting the random ones by tt. Using the 
appropriate indexation of particles (see (1181) ). the calculation of growth rates amounts 
to determine the largest eigenvalue of the N/2 x N/2 matrix L defined by 



V(i,i)e l 



NV 



■ > 



2 J 



COS I 



(27) 



where the particle positions {6^} are distributed according to some fo(9), and where 
the 1/N normalization factor has been voluntarily omitted so that the coefficients' 
distribution do not depend on the number of particles. We used a method based on 
Random Matrix Theory to calculate the growth rate's expectation. When the random 
coefficients verify (Lij) ^ = \i > 0, an extension of Wigner's law [7] states that the 
largest eigenvalue is asymptotically approximated by 



v 



N 



N/2 



+ 



a 



// 



+ o 



(28) 



where a 2 = (L?) . . Since (Ai)j = 1, Theorem 2 in Reference [7] states that v has a 
normal distribution of expectation 1 + (N/2 — 1) \i + a 2 //! and bounded finite variance 
2a 2 . Using Equations dHJ and (J22D, 

one finds the mean squared growth rate as 



' <2 



1 



a 

AH 



(29) 



The value of A 2 is hence distributed according to a normal law of variance 8a 2 /N 2 
because of the 2/N rescaling factor. The expectation of the growth rate reads then 



(7) = AT- 



x cxp 



N 2 
16a 2 



(x - <A*» : 



dx, 



(30) 



where M is the normalization factor. 
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This approach simplifies the calculation of the growth rate, since no more effort 
in the diagonalization of matrix A has to be done. However, the most restrictive 
applicability condition of Equation ( 1301) is the particular symmetry of the initial state, 
which allows the use of Equation ( 122|) . The conditions on the probability distribution 
of Lij are not as limiting as the required symmetry of the initial state. 

4-1. Random uniform bi-cluster 

The equilibrium configuration generated by two symmetric waterbags of parameter AO 
yields a state topologically close to Eqs. ( ITBl . Therefore, we expect the growth rates for 
these random states to be close to the ones given in Equation (1251) . However, the latter 
Equation cannot be the asymptotic form of Equation ( J30l) since the highest eigenvalue 
fluctuates around 1 + (N/2 — 1) fi + a 2 / fi, which corresponds to the eigenvalue of the 
deterministic matrix Oy = fi for i ^ j, an = 1, that is completely different from the 
deterministic matrix L defined by Equation ( 1201) . 

One needs to verify the applicability of Equation fl28|) before calculating the 
expectation of the growth rate with Equation (1301) . fi and a 2 are accessible through 
the following formulae 

fi = (cos (9i - 6 ] )) fQ = jj / (^)/o(^) cos (9i - 9j) d6 l d6 J} (31) 

a 2 + /i 2 = (cos 2 (Oi - 6 3 )) h = jj /o(^)/o(^) cos 2 (9 t - 6 3 ) dft t d%. (32) 

Writing x(N) the characteristic function of the set X, the waterbag distribution of 
parameter A8 has the probability density 

/o(0) = ^x([-A0,A0]). (33) 

Hence, we have 

[i = sine 2 (A0) , (34) 

a 2 = \ + COsH 4 A6) sine 2 (AO) - /A (35) 

Clearly, /i > for A8 < ix j2 and o 2 is finite, which means that Equation (!30l) holds. 
Figure (jlj) shows the numerical fits for randomized bi-clustered initial configurations and 
the theoretical mean value given by Equation fl30|) . We also plotted the growth rate (1251) 
corresponding to the deterministic bicluster to show that, as expected, its behaviour is 
close to the expectation of the random one for a wide range of AO. 

4-2. Random Gaussian bi-cluster 

In this subsection we show another example of use of Equation ( J30l) for a more difficult 
case. The particles are no longer distributed according to a waterbag density, but with 
a Gaussian one. We define fo(0) by 



Me) 



a e erf ( — ) 

\2a e ^2) 



6XP I "2a 2 ' ' (36) 
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Figure 4. Plot of the growth rates from the uniformly randomized bicluster initial 
configuration with respect to A9. The circles correspond to averaged exponential fits 
from the numerically integrated magnetization, the solid curve corresponds to Equation 
(1501) . The dashed curve corresponds to Equation (1251) . Each numerical calculation has 
been performed with 1000 particles, and each circle results from the average of 8 runs. 
We clearly see a very good agreement with the theoretical prediction. 



The normalization factor has been calculated so that the particles are distributed on 
[— 7r/2,7r/2] with a standard deviation Oq. Moreover, 



a 



erf 



1 e 
2 + " 



9^ erf 



7r — 2iat 



erf 



2a e \/2 
-i -l 



> 



(37) 



(38) 



As shown on Figure (jSJ), the agreement between the experimental average and the 
random matrix theory is very good. 




Figure 5. Plot of the growth rates as a function of ae for the random Gaussian 
bicluster equilibria. The circles correspond to averaged exponential fits from the 
numerically integrated magnetization, while the curve corresponds to Equation (|30[) . 
Each numerical calculation has been performed with 1000 particles, and each circle 
results from the average of 8 runs. Here again, the agreement with the random matrix 
theory prediction is very good. 



In the general case, relaxing the symmetry assumption in the preparation of 
the finite-iV equilibria yields a vanishing average value of the matrix elements of A 
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consistently with a vanishing magnetization. This prevents the application of the 
Theorem 2 of [7j and Equation f l30|) no longer holds. To overtake this difficulty, we 
now turn to the more usual continuous approach. 

5. Linear theory within the continuous approach 

Writing f(6,p,t) the distribution function, one has the formal Vlasov equation 

One can define the density of particles n(8, t) and the mean velocity field v(8, t) through 

/+oo j r>+oo 

f(e,p,t)dp and v(6,t) = —— pf{B,p,t)dp.{40) 
n[0,t) 

The HMF force field E(6,t) is then given by 

sin (a - 9) n{a, t)da. (41) 

-7T 

The early evolution of the cold HMF can be reduced to a fluid description. The hierarchy 
of the moments of the Vlasov equation can be stopped at the first order since the 
temperature vanishes. By taking the moment of order zero, one immediately obtains 

Multiplying Equation ( 15§j) by p and integrating over p, one has 

+ I^ p2 % dp ~ mt)nv = °- (43) 

Injecting Equation ( 142]) in Equation fj43]) and dividing by n yields 

dv dv 

_ + „_ -*(,,») = 0, (44) 

in the cold case, for which the mean square of the momentum equals the square of 
the mean velocity. Considering the stationary solution given by some ^o($) yielding a 
zero magnetization, i.e. having a zero m = 1 Fourier component, and v = 0, one puts 
v = 5v(8) exp(iut) and n = n {6) + Sn(8) exp(iwt) in Equations (jUJ), fl42]) and ( 1441) . 
Expanding 5v (8) and Sn(8) in Fourier series, one obtains the linear system 



iuj 5n m exp(im6 l ) + i (m + /) n 0tm 5vi exp [i (m + /) 8] = 0, 

m ml 

\u 2j 5w m exp(im6 l ) + i7r 2J 5n m (5_i im e~ l61 - 5i, m e l61 ) = 0. 



m m 



This is 



w5n fc + k 22 n , m 5v k - m = 0, (45) 

m 

71 

iu5v±i = ± — 6n±i, (46) 
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and 5v m = for m/±l. This gives finally 

u 2 5n k + kn (n ,jfc_i5ni - n 0)k +iSn^i) = 0. (47) 

The dispersion relation is thus given by det M(oj) = 0, where M is generally an infinite 
matrix of elements 

M M = u 2 5 ki i + krc (<^,in ,fc-i - <^,-in 0) fc+i) , (48) 
with (h,£) G Z 2 . 

Let us consider the finite-size 2N+1 square matrix M^ N \ which coefficients coincide 
with M tj , V(z,j) G {-N,Nf. For N= 1, we have 

,2 



" 7™0,0 




-7m _ 2 




-7m 0) 2 w z + 7rn 0i0 

The condition det M = is then fulfilled when 

(w 2 + 7m , ) 2 = n 2 n 0t2 n _ 2 , (50) 

but since n(8,t) is real, no,-2 = ^2 where the bar denotes the complex conjugate. 
Therefore, one can take the square root in Equation (1301) . leading to 

w 2 = -7m ,o ±vr |n 0)2 | • (51) 

It is then easy to show by recurrence that VJV, detM^ 1 ) = det , so that 

Equation (|5T]) gives the only non-vanishing roots to the general dispersion relation. 

This allows to obtain the growth rate for the cold unmagnetized HMF in the infinite 
N limit as 

+ 27r|n 0i2 | 



(49) 



7 



V2 



(52) 



where we used the fact that uq{6) is normalized giving rio,o = 1/27T. 

It is easy to check that the iV — > 00 growth rate (|2^|) may be obtained from 
Equation (|52|) in the uniform bicluster configuration, for which 27rno,2m = sine (2m A0) 
and no,2m+i = 0, Vm e Z. The formula (|52|) was successfully tested for a variety of 
spatially inhomogeneous equilibria with M = (See the plots of the instantaneous 
growth rates M/M with respect to time in Figure [6]). 

We shall now briefly discuss the large-time nonlinear features of the HMF model in 
the light of Lynden-Bell's picture of violent relaxation. 



6. Final discussion on the cold HMF case : an example of violent relaxation 

As two-body, collisional, relaxation timescales were far too large to account for the 
observed luminosity profiles of elliptical galaxies, the concept of violent relaxation was 
introduced by Lynden-Bell in his famous 1967 paper [SII2|- Following this seminal study, 
there appears to be nowadays a common agreement on the scenario of the relaxation 
process for iV-body long-range systems. It is supposed to be divided in two parts: 
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Figure 6. Plot of the numerical growth rate M/M with respect to time for different 
initial equilibria. When the magnetization is almost exponential, M/M is practically 
constant and matches the theoretical value given by Equation (|52l) . represented by 
the horizontal lines. We clearly see that when the second harmonic of no is zero, the 
growth rates are equal to l/y/2 as in the homogeneous case. All the runs have been 
performed using 10 particles. 

the properly speaking violent relaxation part, namely a rapid evolution of the system 
towards a quasi-stationary metastable stable, followed by a much slower thermalization 
phase towards the state predicted by equilibrium statistical mechanics. The HMF model 
provides a simple yet nontrivial long-range system to consider this issue. 

Extensive numerical studies of the HMF model are now available in the literature, 
mostly for two types of initial states: waterbag initial distribution functions in positions 
and momenta and Maxwellian initial distribution functions in momenta with possibly 
waterbag initial distribution functions for positions. Long-time discrepancies between 
time averaged observables and their ensemble predictions have been reported only with 
initial waterbag distributions in momenta, in connection with the emergence of so- 
called quasi-stationary states. Linear theory, mostly done within the Vlasov framework 
and for spatially homogeneous states, can be used as a guideline to discuss relaxation 
properties in the spirit of Lynden-Bell's picture. For instance, the spatially homogeneous 
cases with initial waterbag distributions in momenta are unstable with growth rates 
equal to a/ 1 /2 — 3T where T is the initial temperature associated to the waterbag. 
Consequently, for T > 1/6 the system becomes linearly stable. This case corresponds to 
a density energy U equal to 7/12 [3j[T0]. This value happens to coincide with the energy 
threshold value above which pathological relaxation behaviours have been reported for 
spatially homogeneous, waterbag in momenta, initial distributions functions. There 
is then a simple, but apparently not explicitly formulated, manner to interpret the 
relaxation features observed in the HMF model: an efficient, rapid, convergence towards 
equilibrium statistical mechanics' predictions requires starting from the vicinity of an 
initially unstable equilibrium. Linear instability then naturally provides the violent 
relaxation phase. The cold case, on which we have just focused, illustrates well this 
scenario. Actually, large-time simulations show that the modulus of the magnetization, 
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for instance, quickly converges towards a saturated state with a value that is close to 
its ensemble prediction. It is interesting to note that similar observations have been 
reported for another long-range, mean-field, wave-particle Hamiltonian model starting 
with a cold beam of particles [12J. An heuristic way to explain this rapid convergence 
is through the comparison between the linear and nonlinear timescales available in the 
system. The linear timescale is obviously given by the e-folding time, namely by the 
inverse of the linear growth rate 7. The nonlinear timescale is the trapping time of the 
particles in the mean-field potential well, namely the inverse of the bounce frequency 
ojb = M 1 / 2 . Nonlinear saturation takes place when both timescales balance, namely for 
7 ~ Ub- Therefore, one expects nonlinear saturation to take place close to equilibrium 
predictions provided that the linear growth rate is of the order of the ensemble average 
of the nonlinear frequency. 

This is indeed the case for the cold system for which U = 1/2. According to 
equilibrium statistical dynamics [IT] , this gives (M) M = (M) c = l/\//3 = vT where the 
inverse of the temperature is given implicitly by Ii(a/5)/Io(a/5) = l/v% Numerically, 
the ensemble average of the magnetization for the cold HMF is then (M) c ~ 0.62 and 
the equilibrium temperature is T ~ 0.39, which does yield 7 ~ (w&) M . Figure [7] shows 
that the cold unmagnetized HMF does experience a violent relaxation. The figure makes 
however apparent that that this phase is followed by a much slower thermalization phase 
that is needed for phase space sweeping and convergence towards equilibrium statistical 
predictions. Actually, in the early nonlinear saturation, the magnetization oscillates 
about 0.60, below the ensemble equilibrium value. The drift towards the equilibrium 
statistical predictions takes place on a much longer timescale than the violent relaxation 
timescale, which is the nonlinear saturation timescale of order 7 -1 . 



0.8r 




12 3 4 



iog 10 (t) 

Figure 7. Time evolution of the magnetization starting from a gently perturbed 
homogeneous quiet start of the cold HMF. 
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